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Abstract. In the present study we examine non-Gaussian spreading of solutes 
subject to advection, dispersion and kinetic sorption (adsorption/desorption). 
We start considering the behavior of a single particle and apply a random walk 
to describe advection/dispersion plus a Markov chain to describe kinetic sorp- 
tion. We show in a rigorous way that this model leads to a set of differential 
equations. For this combination of stochastic processes such a derivation is 
new. Then, to illustrate the mechanism that leads to non-Gaussian spreading 
we analyze this set of equations at first leaving out the Gaussian dispersion 
term (microdispersion). The set of equations now transforms to the telegra- 
pher's equation. Characteristic for this system is a longitudinal spreading, 
that becomes Gaussian only in the long-time limit. We refer to this as ki- 
netics induced spreading. When the microdispersion process is included back 
again, the characteristics of the telegraph equations are still present. Now 
two spreading phenomena are active, the Gaussian microdispersive spread- 
ing plus the kinetics induced non-Gaussian spreading. In the long run the 
latter becomes Gaussian as well. Another non-Gaussian feature shows itself 
in the 2D situation. Here, the lateral spread and the longitudinal displace- 
ment are no longer independent, as should be the case for a 2D Gaussian 
spreading process. In a displacing plume this interdependence is displayed as 
a 'tailing' effect. We also analyze marginal and conditional moments, which 
confirm this result. With respect to effective properties (velocity and disper- 
sion) we conclude that effective parameters can be defined properly only for 
large times (asymptotic times). In the two-dimensional case it appears that 
the transverse spreading depends on the longitudinal coordinate. This results 
in 'cigar-shaped' contours. 

Keywords;Advection-diffusion equation and kinetic adsorption and random 
walk and Markov chain and solute transport and telegraph equation. 



1. Introduction 

It is well known that in the field a contaminant plume spreads at a higher rate 
than as predicted by theory and lab experiments. In addition, one observes that 
the spreading pattern often deviates from the Gaussian pattern, especially at early 
times (tailing). These phenomena have been analyzed intensively, theoretically, 
numerically and by field experiments ([in],[Mj,I7],[II],[Il]),[3D],[3],[ll]). The 
majority of papers on this subject attribute these phenomena to the heterogeneity 
of the medium. In our paper we show that in a homogeneous aquifer a similar 
behavior occurs when contaminants are subject to a relatively slow kinetic adsorp- 
tion/desorption reaction. 
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We consider a homogeneous medium and follow a (solute) particle during its 
motion through the pore system, while simultaneously it is subject to adsorp- 
tion/desorption. Although the porous medium itself is simply homogeneous and 
non-stochastic, the particle's behavior is still chaotic and, therefore, our particle 
model will be stochastic. Note that our approach differs from that of papers on 
particle tracking in the sense that we do not start with the advection-dispersion- 
reaction equations and, accordingly, perform a particle tracking algorithm to solve 
or simulate the equations. Instead, we start with a stochastic description for the 
movement of a single particle and work our way towards the differential equations. 
Particle tracking papers assume an analogy between particle tracking and the dif- 
ferential equations either by just stating the validity of this analogy or by referring 
to previous papers on particle tracking. To our knowledge a rigorous derivation of 
the differential equations for advective-dispersive transport with kinetic sorption 
starting from a stochastic model for a single particle does not exist in any of those 
papers. 

We also discuss particle models used in various other fields, such as the velocity 
jump model (chemotaxis or movement of bacteria, Hillel). This velocity jump 
model (Giddings-Eyring model) leads to a telegraph equation. It differs from our 
model by the fact that there is no diffusion or dispersion phenomenon. We will 
demonstrate that, when the velocity jump is extended with dispersion/diffusion, 
the telegraph character stills remains. This explains the non-Gaussian behavior in 
the pre-asymptotic stage. 

The non-Gaussian features can also be illustrated by (spatial) moments, espe- 
cially by the third (skewness) and fourth (kurtosis) centralized moments. For the 
ID situation first and second moments were studied by previous authors e.g. |34j . 
[5] and |31| . These papers mainly focused on the large time (asymptotic) results. 
Pre-asymptotic expressions for first and second moments were first derived by |27j . 
Unfortunately, in the formulation of the solutions Michalak and Kitanidis suggest 
that for an arbitrary initial distribution of the phases the mean and variance in each 
separate phase can be obtained by a linear combination. However, this is only true 
for the non-centralized moments, but not for the variance. Moreover, one of the 
4 solutions is incorrect (see section 5.1, also [lO]). In our paper we shall consider 
the moments for the 2D situation and examine in detail the conditional moments. 
These moments manifest more clearly the non-Gaussian behavior than the ID or 
marginal moments. It shows that non-Gaussian behavior can exist, even when the 
marginal moments suggest that the behavior is Gaussian. 

In our 'particle view' the behavior of a particle is described by a random walk 
and a two state Markov process. A similar model was used by [34], [35] and [37] . 
Equivalent processes occur in other fields, such as chromatography ([H], [TB], [35]) 
or in statistical physics (e.g. [13], [TS], [H], [33] and [2B]). The amount of literature 
related to this topic is extensive and spreads over many different fields. Here, we 
look at this topic from the perspective of solute transport in groundwater. 



2. The particle view 

2.1. A stochastic particle model. We model the movement of a single particle 
subject to an advection/dispersion/sorption process over a time interval [0,i]. For 
simplicity we discuss the one-dimensional case. We discretize time by choosing 
some n, and by dividing [0,t] into n intervals of length At — t/n. We observe the 
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state of the particle at the tnne pomts 0, Ai, 2At, . . . , nAt. For this state there are 
two possibihties: 'free' or 'adsorbed,' which we code by the letters / and a. The 
particle can only move when it is 'free,' and its displacement has two components: 
dispersion and advection. 

Let Xk be the displacement due to the dispersion of the particle the kth time 
that it is 'free.' We model the as independent random variables with mean and 
variance 

E [Xk] = and Var (Xk) = 2D At. 

The displacement due to advection is given by vAt, where v is the (deterministic) 
advection velocity. 

Let Kn be the number of intervals [kAt, [k + 1) At) during [0, t] that the particle 
was 'free.' (Here the open bracket at the right end indicates that the point (fc + l)At 
is not included.) In other words, KnAt is the free residence time of the particle in 
[0,t]. Let S{t) be the position of the particle at time t. Combining the two types 
of displacement we obtain 



fc=i 



vAt). 



The distribution of Kn is determined by the kinetics, i.e., by the switching between 
the 'free' and the 'adsorbed' state. This is naturally described by a two state 
Markov chain. The state-transitions of this chain after a certain time step At are 
given by a transition probability matrix (pij): 



' Pff 


Pfa ' 




Paf 


Paa 





1 



1 



This means that for instance the transition from 'adsorbed' to 'free' has probability 
a = Paf to happen during the time interval [kAt, {k + I) At) (note that actually we 
make this change — if it takes place — at the end of the interval) . 

We will be mainly interested in the moments of S(t). Below we will compute 
the first and second moment, and in the next Section we discuss the (centered) 
third and fourth moment. To compute E [S{t)] we use the well known formula (see 
e.g. [29 ) for a random sum of Kn independent and identically distributed random 
variables Yk (also independent of Kn): 



E 



'K„ 

E 

.k=l 



Yk 



E [Kn] E [Y{ 



Here the mean of Kn equals: E [Kn] 



This expression can be obtained 



vt. 



from [31], or [TU]. Substituting we find (with nAt — t) 

E [Sit)] = E [Kn] (E [Xi + vAt]) = - 

a ^ u 

Here, (a + b) /a is the retardation factor R. To see this, note that the probability 
vector (b/{a + b) a/{a + b)) is the stationary distribution of the Markov chain, and 
so a/{a + b) is the fraction of time the particle is free. Thus the effective velocity 
V* = v/R. We have implicitly required that the particle at time is given the state 
'adsorbed' or 'free' according to this distribution, for other initial distributions there 
will be a correction term in the formula for E [S'(t)], which tends to as t tends to 
infinity (cf. [TU]). 
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To compute Var(S'(<)) we use the well known formula (see e.g. |29]) for the 
variance of a random sum of i.i.d. random variables Y^. (also independent of 



(1) 



Var ^ Ffc = E [X„] Var (Fi) + Var (/f„) (E [^1])^ 



This yields with Yk = Xk + f Ai and nAt = t: 



(2) 



Var {S{t)) ^ E [K,,] Var {Xk + uAt) + Var (X„) (E [Xfe + vM]f 
2Dt + Ysir{Kn)v^{Atf, 



a + b 



where (as can be deduced from ^S] or [TU] ) 



(3) 



Var(^C„) 



ab{2 -a-b) 2ab{l - a - b) 



{a + by 



{a + 6)4 



[l-{l-a-br] 



The equations ^ and ([3| thus tell us that the variance of the displacement of the 
particle grows more or less linearly in time with (asymptotic) slope 



o „ o,b(2 — a — b) 9 . 

2D H \ rri^V^At. 



a + b 



(a + by 



2.2. Skewness and kurtosis. To obtain the skewness of S{t) we must use the not 
so well known formula for the third central moment of a random sum of Kn i.i.d. 
random variables Yfc (also independent of Kn). 



E 



{S{t)-E[Sit)]] 



= E 



K„ 



fc=l 



.fc=i 



= E [Kn] E [{Yi - E [Yi]y] + 3E [Fi] Var (Fi) Var (if„) 
+ (E [Y,]yE[{Kn-E[Kn])'] . 

To actually derive a formula for the skewness of the displacement of the particle 
from this equation will lead to very heavy computations (and the situation for 
the kurtosis is even worse). However, without doing any computations we can 
already tell that as t — >■ 00 the skewness must tend to zero, and the kurtosis to 3: 
this is because the distribution of the displacement of the particle will tend to a 
Gaussian by the Central Limit Theorem for random sums of independent identically 
distributed random variables. In our case this follows since A'„/n tends in the mean, 
and hence in probability to a/ {a + 6), see [T^, page 258. 



2.3. Decreasing the time steps. The discrete time steps are somewhat unnat- 
ural. We would like to let At tend to 0. But then we have to realize that a and 
b are functions of At. Since the probability that the particle changes its state is 
proportional to the time At it is observed (if At is not too large), we should put 



a = liAt, b = XAt, 
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where /i and A are now the rates at which the particle switches from 'adsorbed' to 
'free', and from 'free' to 'adsorbed'. Substituting this in Eqs ^ and ([s]), we obtain 



Var(S'(i)) 



-2Dt 



Am(2 - (A + Ai)AO 



A + ^ " ' (A + //)3 
2AAt(l - (A + /i)At) 



Letting At — 0, and hence n — >■ cx) we obtain 



1 



l-(A + /i)- 



Yar (S(t)) ^ —^2Dt 



2A/i 



2A/i 2, 

(A + mP" (A + ^)4 



-(A+M)t 



Thus we recuperate a (more general and more detailed) version of the main result 
of and there is a match with the expression that comes from the moment 
analysis based on the differential equations (as can be derived by correcting the 
results in [27]). 



2.4. The state of the particle at time t. Let Sf{t) be the position of the particle 
at time t given that it is free at time t. 

To find the distribution of Sf{t)^ we need the distribution of Kn \ the number of 
intervals [fcAt, (fc + l)At) during [0, t\ that the particle was free, given that it is free 
at time t = n/S.t. We find now (where E 



E[5/(t)]=E 



vAt 



can be deduced from |35]) that: 
6(1 - (1 - a - 6)")" 



a + b (a + 6)2 

b(l-(l-a-6)*/^*) 

, ' 7 TT^ vAt. 

a + b (a + 6)2 

Substituting a — jiAt, b = AAt, and letting At — !• we obtain 



vAt 



vt- 



vt - 



Xv 



A+M (A + Ai)2 
From |10[ we have that Var (^Ki^^^ equals 
"a6(2-a-6) 2 6 (a - 6) (1 - a - 6)' 



1-e 



-(A+M)t 



(a + 6)3 (a + 6)3 

"6(3a-6) 4a6 



(a + 6)3 (a + 6)4 
Using Equation ([ij we derive from this 



[1 - (1 - a - 6)"] + 



(a + by 



[l-{l-a~bf 



Var(5y-(t))=E 



if) 



Var {Xk + vAt) + Var ( K^/^ ) (E [Xk + vAt]) 



a + b 



2Dt + Var 



(i^(/)).2(At)2 



6(l-(l-a-6))*/^* 
(a + 6)2 



2i:)At. 
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Substituting a — fJ-At, b — XAt, and letting At — > we obtain 



Yar {Sf{t)) = —^2Dt 
A + /i 

2A^ 



A 

2A(^^-A) (^+^), 



(a + m)3 ix + ^ir 



4A/i 



1 



,-(A+M)t 



A2 



1 - e 



It can be shown that this matches with the expressions in |27j . when these are 
corrected as in |10| . Similar computations can be made for the displacement of the 
particle given that it is absorbed at time t. 

2.5. Derivation of the diflFerential equations. We will now show how the fun- 
damental differential equations (10 1 can be obtained from a diffusion limit of 
the single particle model. Our approach is similar to the one followed for transport 
in fluidized beds in [3]. In order to obtain this diffusion limit we also discretize 
space in locations 

iAx, i = . . . , —1, 0, 1, ... . 

Here we let Ax depend on At in the classical way, which is motivated by the fact 
that typically at time t the spatial fluctuations arc of order 



(4) 



Ax = c\/At, 



where c > will be chosen later. The particle moves according to a Markov chain 
(Zn), which is a birth-death process (birth=one step to the right, death=one step 
to the left), with the additional possibility that the particle may become adsorbed 
and free again. The state space is therefore a product 

5-{...,-l,0,l,...}x{a,/}, 

where e.g. Zn ~ {i, a) means that at time nAt the particle is at iAx, and is 
absorbed. Let 6i and be the probabilities that the particle (in the free state) 



5r 




{i - l)Ax iAx {i + l)Aa 



Figure 1. The three possible movements of the particle when it is free. 

moves from iAx to {i + 1) Ax, from iAx to {i — 1) Ax, or stays at iAx. Note that in 
this model the dispersion may depend on the location, i.e., we could more generally 
consider (deterministic) inhomogeneous media. To make this more realistic we 
could also let the adsorption and desorption rates depend on z. As we consider here 
only the homogeneous case, we take j3i = /3,Si — S and Ui — a for each i. To fit 
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the new model with the model in Section [2. 1| we choose values for S and a in the 
following fashion: 

At v"^ / At\'^ V At 
' XAx^ ^ ^\Ai) ^ 2 

At w2 / At \2 V At 
JAx^ ^ ~2\Ai) ^ 2 

(Aa;)2 \AxJ 



By fitting we mean that the mean and the variance of the displacement in a time 
interval of length At of the particle (in the free state) are the same in the two 
models. Indeed, the mean of this displacement equals /3Ax — SAx — vAt, and the 
variance equals 

PiAx)"^ + SiAx)"^ - {vAtf = 2D At. 
In terms of At only, using Q, the displacement probabilities are 

D v^At vVAi 



^ " c2 ' 2c2 ' 2c 
D v^At v^/At 

_^ 2D v'^At 



From this we see that these are indeed probabilities for At small enough, provided 
we choose c > \J2D. The possible transitions of the chain are 




(2-1,/) ^ (i,/) ^ + 



(i — 1, a) ,. , (j + 1, a) 




The corresponding transition probabilities are: 

= - AAt) P(j,/),(i+i,a) = /?AAt P(tj),{i-ij) = - AAt) 

P{i,f).ii-l.a) = (^AAt P{t.a).ii,a) = 1 " t^At P{t,a),{i.f) = A«At 

P{iJ)AiJ) = Q^Cl ^ -^At) P{i.f),{i,a) = aAAt 

Let pf{i,n), respectively Pa{i,n), be the probability that at time t = nAt the 
particle is free, respectively absorbed and at position lAx. 
The master equations for the particle are 

Pf{i,n+1) =p/(i- l,n)/3(l - AAt) + + 1, n)<5(l - AAt) 
+ Pf{i,n)a{l — AAt) + pa{i,n)fiAt, 
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and 
(6) 



Paii,n+ 1) ^pf{i- l,n)l3XAt+pf{i + l,n)6XAt 
+ Pf{i, n)aXAt + Pa{i, 'n){l — fiAt). 



With some abuse in notation regarding the functions Pa and p/ we wiU denote 
the limiting probabihties of Pa{i,n) and pf{i,n) as lAx x, and t = nAt by 
Paix, t) and Pf{x, t). We wiU obtain partial differential equations for these limiting 
probabilities when we let At 0, and i = i{At) oo, in such a way that iAx — ?► x. 
(The obvious way to achieve this is to take i{At) equal to the integer closest to 
x/{c\fAt).) Rearranging ([s]) we obtain 

Pf[i,n + I) — pf{i ~ I, n)j3 + 1, n)5 + Pf[i, n)a — \At[pf{i — 1, n)j3 

+ Pf{i + l,n)6 +pf{i,n)a] + pa{i,n)fiAt. 

In the first three terms of the right hand side we substitute the values for /3,6 and 
a: 

(8) 



Pf[i 



l,n)/3 + pf{i + l,n)S + pf{i,n)a 
At w2/Ai\2- 



D 



(Ax)2 
V At 
2Ax 



2 \AxJ 



Pf{i - l,n) - 2pf{i,n) +Pf{i + l,n) 



Pf{i~l,n)+pf{i + l,n) +pf{i,n) 



Substituting Equation ([8| into ([t]) and dividing by At, we obtain an equation for 
the difference quotient {pf{i + l, n)—pf{i, n))/ At. Then, letting At 0, we obtain 



(9) 



dpf{x,t) _ ^d'^pf{x,t) _^ dpf{x,t) 



Xpf{x,t) + npa{x,t). 



dt dx'^ dx 

For the adsorbed phase a similar equation can be derived. Here, since the adsorbed 
particle is not subject to advection or dispersion, the first and second term at the 
right hand side are absent. 

dpa{x,t) 



(10) 



dt 



'fiPa{x,t) + Xpf{x,t). 



2.6. Prom particle to plume. It might seem surprising that we study the behav- 
ior of a contaminant plume from the stochastic analysis of a single particle. Here 
we illustrate how these two approaches are connected. We model the contaminant 
plume by a collection of N particles. These particles move independently according 
to the same law as the single particle considered in the previous sections. Let Si{t) 
be the position of the ith particle at time t. We are interested in the centroid Z{t) 
of the plume at time t. This is given by 



Z{t) 



1 ^ 



i=l 



We are also interested in the spreading of the plume around its centroid. This we 
measure by the (empirical) variance V{t) of the particles given by 



(11) 



1 ^ 
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The random variable Z{t) is an average of independent identically distributed 
random variables with finite expectation. Therefore by the strong law of large 
numbers for t fixed, and N large 



(12) 



Z{t)^E[S{t)] 



We now turn to the spread of the plume. Here the situation is more complicated 
because the terms in the sum are no longer independent random variables, and also 



depend on N. It is well known that a rewriting of Equation (11 1 yields 
1 ^ 



V{t) 



We have already argued (see ( 12 )) that the last term is approximately 0, and another 
application of the law of large numbers to the first term yields that 

nt)«E[(5(0-E[5(t)]f] =Var(5(i)). 



3. Giddings-Eyring Model 



A simpler system, different from but still related to the system described by Eqs. 
^ and ([To]), is given by the following equations: 



dNf{x) dNf{x) , , ^ N 

= XNf{x)-^^N^{x). 

The difference is that the dispersion process is absent. This system is known in the 
literature under various names. In probability theory and statistical physics it is 
known as a persistent or correlated random walk and has been studied e.g. by |14j . 
[H]) [33], HB) and [32] ■ In the field of chromatography the system is intensively 
studied as well ([T7], [T5], [Hj). Solutions for the probability density functions of 
the particles are given by Giddings and Eyring. [50], citeHillel applies the equa- 
tions to the movement of bacteria in the direction of the gradient of food molecules 
(chemotaxis) and uses the term "velocity jump process" . Several interesting obser- 
vations can be made with respect to the spreading of the particles. The particles 
are undergoing an 'apparent' dispersion, despite the fact that (hydrodynamic) dis- 
persion is not included in the model. This 'kinetics-induced' dispersion develops in 
a non-Gaussian way. In systems where both hydrodynamic dispersion and kinetics- 
induced dispersion are present, the latter sometimes can be more dominant, such 
that the non-Gaussian spreading is observed also in systems with hydrodynamic 
dispersion. We shall illustrate this later. First, we introduce a moving coordinate 
system with velocity v* . The new x-coordinate is: 

X = X — v*t, 

where v* = vfj,/(X + /i). In the new coordinate system free particles move to the 
right with velocity Vf ^ v — v* ^ vX/(X + fi). The adsorbed particles 'move' with 
velocity Va — —v* — —v x fi/{X + The minus sign indicates that the movement 
is to the left (i.e. with respect to the new coordinate system). The equations are 
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now: 



dNf{x) dNf{x) , , ..r,^ 

dNaix) dNgjx) 

+ Va = >^Nf(x) - ^lNa[X) 



dt 



dx 



The equations can be rewritten as a single differential equation by considering first 
the sum and difference of the particle distributions (Kac's trick), i.e., 

u{x) = Nf{x) + Na{x); w{x) = Nf{x) - Na{x). 

After summation and substraction of the differential equations we obtain for u and 



du 
dt 
dw 

Ik 



X — fi\ V du V dw 
X + fi) 2d^ ~^ 2~d¥ 
X — fi\ V dw V du 
X + fi) 2~d¥ ^ 2dS; 



= 



'u{X — /i) — w{X + /i). 



Now we differentiate the first equation to t, the second to x and eliminate the 
derivatives of w: 



(13) 



Xfiv"^ d'^u 1 d'^u V (A 



(A + fif dx^ A + ^ dt^ (A + p,f dtdx 



ji) d^u du 

" Ik 



This is a telegrapher's equation with an additional term due to asymmetry {fi ^ X). 
Hillel discusses the symmetrical (A = fi) case only, but Weiss ([SI]) also makes some 
remarks on asymmetry. Also see Masoliver [26] and [5]. The telegraph equation 
may be interpreted either as a diffusion equation with a perturbation term that 
disappears at large times, or as a wave equation with a perturbation term that 
disappears at early times. Thus, as time proceeds, the system can be described by 
three different equations: first, a wave equation for early times; secondly, a telegraph 
equation for intermediate times and thirdly, an advection-dispersion equation for 
large times. 

3.1. Large time behavior and the advection-dispersion equation. Hillel 
uses the following argument that leads to a useful result for large times. For large 
times the velocity [LT~-^] and sorption rates [T~^] typically are expressed in large 
time units and thus their values become large. Then, the first term at the left 
hand side of ( 13 ) dominates over the second and third term. Therefore, at large 



times the equation approximately describes a dispersion process with an equivalent 
dispersion coefficient D* , purely induced by the kinetics: 



D* 



A/iw^ 



(x+f^r 



3.2. Short time behavior and the wave equation. In a similar way it can be 
shown that for small times the terms at the left hand remain, while the right side 
becomes small. The remaining expression is a wave equation that can be written 
as: 



'd^ 
dt 



fiv 



d 



X + fj. dx 



' d^ 
dt 



Xv d ' 

X + ^ dx 
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For an initial pulse at x = the solution represents two pulses propagating along 
the characteristics. In x-space, 

X — Vat = 0, X ~ V ft — 

which in x-space correspond to: x = 0, x — vt = 0. 

These pulses can be identified as a stagnant pulse (the original adsorbed parti- 
cles) and a travelling pulse (original free particles). 

3.3. Intermediate time behavior. The wave equation and diffusion equation are 
approximations for the process at short and large times respectively. The interme- 
diate time is described exactly by the full telegraphers's equation. Therefore, exam- 
ination of the solutions to this equation will give information on the pre-asymptotic 
spreading behavior of this process. 

Solutions for the distribution of the particles have been derived by [17], [16], 
and [22, . For detailed discussions of these functions see [36] , [23] and [37] . 

Giddings and Eyring consider four types of densities hff, haf, hfa, haa'- 



hffir, t) = e-^--^(*-^) ^ ^/i {0) + e-^'5{t - r) 

hfa{T,t)^\e-^--^^'--ha{e) 

haf{T,t)=f^e-^^-^^'--hoi9) 



(14) 



h^{r,t) = c-^r-,(t^r)^ Mt^ r) ^^ {e)+e-^'5{T), 

where r = x/u, 9 = 2-\/A/ir(i — r) and and /i(-) are modified Bessel functions. 
Note that t = x/v \s not simply a convenient scaling of the x-coordinate, but t 
also represents the residence time in the free phase. The expressions hij represent 
the probability densities of the free residence time for different phases and different 
initial states of the particles. The first index indicates the initial state of the 
particle and the second index indicates the state of the particles the pdf is referring 
to. The distributions are zero for r < and r > i (x < and x > vt). The delta 
functions at r = and r = (which in x-space corresponds to x = and x — vt) 
represent exponentially decreasing pulses and can also be identified as the fractions 
of particles that, since t = 0, did not (yet) perform a change of state. 

In Figure [2] and [3] we present graphs with the evolution in time of these distri- 
butions, using as initial condition the (unit) pulse consisting of free and adsorbed 
particles in equilibrium. Let tt/ and iTa be the initial amount of particles in each 
phase. Equilibrium exists for ^ = j, and, if the total amount is unity we have: 



Thus, for t = 0: 



Nf{x,t) = 7r/(5(x); Na{x,t) = 7ra(5(x). 
dis 

hf{T,t) = Trfhff{T,t) +TTahaf{T,t) 
hl^{T,t) = TTfhfa{T,t) +TTahaa{T,t), 

while for the total amount of particles we have, h^gf. = + h'^J^. 



We denote the residence distributions due to this initial condition as /i^' and /i^^: 
(15) 

h7{T,t) = TTfhfaiT.t) +Trahaa{T,t), 
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Figure 2. Distributions of separate phases and total plume at 4 
moments in time. Symmetric case (A = fj,). Solid (blue) line for 
total plume. Dashed lines for (green) adsorbed phase and (red) 
free phase, yellow background for Gaussian distribution. 



and /ifgj with a Gaussian 
distribution with mean v*t and variance 2D*t representing the distribution of a 
solute with velocity v* and dispersion D* . All distributions are plotted versus a 
scaled variable x: 

X — v*t 
" V2D't' 

Several stages of the system are shown. At early times {t — 1/4, Figure [2j upper 
left graph) the pulses of free and adsorbed mass are still distinguishable. The 
pulses move apart in x space and are damped. The mass that 'leaves' the pulses 
gradually fills the space in between and builds up a distribution that becomes 
Gaussian in the end (e.g., t = 16, Figure [2] lower right graph). The distribution in 
the interval between the pulses is absent in a pure wave system and is typical for 
the telegraph equation. In the final stage the pulses are completely damped and the 
distribution approaches the normal distribution. Summarized, at early times the 
'wave-character' dominates, at large times the 'diffusion-character' dominates, while 
at intermediate times the system is adequately described by a telegraph equation 
(travelling and dampened pulses -|- mixed zone in between). For A = /x the Gaussian 
distribution is reached slightly faster than in the asymmetric case (A ^ /Lt), as seen 
in Figure [a] For both cases the pulses disappear for t exceeding both 3/A and 3//x. 
In the next section local dispersion is included. We show that the two-dimensional 
distribution still may deviate from two-dimensional Gaussian functions even though 
the corresponding one-dimensional distribution is close to Gaussian. 



In Figures m and § we compare distributions i 
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Figure 3. Distributions of separate phases and total plume at 4 
moments in time. Asymmetric case (A 7^ /i). Solid (blue) line for 
total plume. Dashed lines for (green) adsorbed phase and (red) 
free phase, yellow background for Gaussian distribution. 



4. GiDDINGS-EyRING model INCLUDING DISPERSION 

We extend the Gidding-Eyring model by including longitudinal and transverse 
dispersion. This way we obtain the following two-dimensional system. 



(16) 



dNf d^Nf d^Nf dNf 
at ox^ oy^ ox 



dt 



\Nf-^iNa, 



where now Nf and Nf represent a two dimensional particle density distribution 
[L~^] and Dl and Dt are the longitudinal and transverse dispersion coefficients. 
We do not use different notations for the one- and two dimensional particle densities. 
We assume that their distinction will be clear from the context. 

It appears that the effect of transverse dispersion is much more dramatic than 
that of longitudinal dispersion, while its solution is much easier to derive. Therefore, 
we analyze longitudinal and transverse dispersion separately and start by including 
transverse dispersion first. 

4.1. Transverse dispersion. We use an approach proposed by [37]. Consider two 
distinct species, one adsorbing and one non-adsorbing. Let the spatial distribution 
of the non-adsorbing solute be given by c{x, y, t). Further, let r be the 'free residence 
time' of the adsorbing particles and let the distribution of r at time t be hij{T,t), 



i.e. for particles in phase j with an initial unit pulse in phase i (see by ( 14 )). If the 



initial pulse is iV°, the particle fraction with free residence time r at time t becomes 
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N^hij{T,t), where t < t. The spatial distribution of this fraction is equal to that 
of the non-adsorbing particles at t — t, or: dNij{x, y, t) — N^hij(T, t)c{x, y, r)dr. 
Summing fractions with r, < t < t, we obtain: 

(17) (x, y, t)^N^ [ (r, t)c{x, y, t) dr. 

Jo 

We apply Van Kooten's approach first to the case with only transverse dispersion. 
For a non-adsorbing solute with advection in the x-direction and dispersion in the 
2/-direction {Dl is assumed zero) the distribution is: 

1 f y"^ \ 

cix, y, t) = — , exp — S(x — vt). 



Insert this function in the integral (17). Because of the delta function the integral 



can be evaluated directly. After substitution of t hy x/v we obtain: 



(18) ''^^^^^y^'^ = '''^M^-^-i-w^. 



For a given value of x, expression (18 1 describes the distribution in y-direction of a 
certain amount of particles. The amount is equal to Nfhij[^, t) per unit of length in 
x-direction and it spreads in the y-direction as a Gaussian distribution with variance 
2Dt-- The transverse variance now depends on the x-coordinate, which is clearly 
in conflict with a 2D Gaussian distribution. The interdependence of transverse 
variance and x-coordinate can be understood by considering the residence times in 
the free phase. During the 'free phase' time particles travel in positive x-direction 
and simultaneously spread in the y-direction. At a given time t the particles that 
have spent more time in the free phase are found further along the x-direction. 
They are also more widely spread in y-direction, since they have been subject to 
dispersion for a longer time. Further, note that for f > i the functions hij{^,t) are 
zero. Therefore, Nij{x,y,t) is zero for x > vt. 

In Figure |4] this is illustrated by 2D contours for unit initial pulses in the free 
and adsorbed phase. The following scaled coordinates are used: 

X — v*t „ y 



(19) 



For a Gaussian distribution with dispersion coefficients D* and Dt/R in the longi- 
tudinal and transverse direction one would expect elliptic contour-lines. Figure |4j 
however, shows typical 'cigar'-shaped contours. It is clear that close to the original 
injection point (in Figure |4] around x = —2) the particles have spend a relatively 
short time, r, in the 'free' phase. During this time the particles displace not only 
very little in longitudinal direction, but also the spread in transverse direction is 
very limited. The result of this is that x displacement and y-spread are no longer 
independent, which is in contradiction to a truly Gaussian system. Also note that in 
the ID case in Figure [2] the same parameters have been used and the corresponding 
distribution in x-direction comes out almost as a Gaussian curve (Figure [2] upper 
right graph). 



4.2. Longitudinal and transverse dispersion. Analytical solutions including 
longitudinal and transverse dispersion have been obtained by [1], but not in a 
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Initially adsorbed 




Initially adsorbed 





Figure 4. Contour-lines for free phase (red) and adsorbed phase 
(green) for several initial conditions. The symmetric case (A = 5 
and /Lt = 5) at time t =1. 



closed-form. The integral expression by Carnahan is, in essence, equal to the one 
obtained by Van Kootcn's approach: 



(20) 



Nij{x,y,t) = 



VDlDtJo 47rr 



Hj(T,i)^_(2 



We apply expression (20) to evaluate the 'full' system (i.e, including advection, 
dispersion and sorption) using the following values: A = 0.2, /i = 0, u = 1, Dl 
= 0.5 and Dt = 0.1. The kinetics induced dispersion coefficient D* now becomes 
0.62, which is in the same order of magnitude as Dj^. 

The contour-lines in Figure [5] represent the distribution of the free phase for t 
— 1, t = 20 and t = 150, for the case of an initial pulse in equilibrium. Note 
that for longitudinal and transverse coordinate we have applied the same scaling 
as in Figure 4, (see (19)). At a short time {t = 1) almost circular contours occur, 
which suggests Gaussianity. At this stage microdispersion dominates the spreading 
process, which now progresses in a Gaussian way. At an intermediate time (t = 
20) the 'cigar'-shaped contours start to develop. Here, we observe an increasing 
influence of the non Gaussian kinetics-induced dispersion. Finally, at large times 
(t = 150) the distribution becomes Gaussian again, with elliptic contours. The 
parameters of the early and late Gaussian distribution are quite different. For short 
times we have velocity v and dispersion coefficients and Dt- For large times 
the velocity becomes v* and dispersion coefficients D^/R + D* and Dt/R for the 
longitudinal and transverse direction respectively. The early Gaussian distribution 
occurs because directly after the start of the pulse the effect of kinetic exchange 
between the phases is still small and the free phase consists mainly of particles that 
did not yet change their state. Therefore, they behave as a non-adsorbing solute. 
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Figure 5. Full two-dimensional case. Distribution of the free 
phase at various times for an initial pulse in equilibrium; sym- 
metric A = 0.2 and ^l = 0.2; Dl = 0.5; Dt = 0.1; v = 1 



As time goes on the influence of kinetics becomes more apparent and 'cigar' shaped 
contours develop. 

5. Moment Analysis 

5.1. One-dimensional. The method of moments was introduced by [2], while 
studying the flow of solutes through a capillary tube (Taylor dispersion) . Later, the 
method was successfully applied to solute transport in a layered porous medium 
by |25) . For contaminant transport with sorption the method was applied by |34j . 
|35) . |S] and [3T]. These authors examine the asymptotic values of zero, first and 
second moments. Solutions for the entire pre-asymptotic regime of zeroth, first 
and second moments are derived by [57]. However, one of their expressions for the 
central second moment turns out to be incorrect. We copy here the formula (the ® 
is a but should be a — ) from [57], page 2136 for the normalized central moment, 
where the solute is in the free phase at time and at time t : 



\ k[l + (iAf{p + lf ) \{p + l){l + pA) 
2v^l3 (1 ~ A){il3^A-i- p{A® 1)) 4D /3 (1 - A) 



k-^{l + l3Af [p+lf k{l + l3A){fi + lf' 

Here Michalak and Kitanidis abbreviate A = A{t) ~ cxp(— (/3 + l)kt), and use the 
notation 

k — ^, /3 = A//i. 

By giving the expressions for aj^, aj^, a^j and a^^ Michalak and Kitanides sug- 
gest that these expressions may be used to obtain variances for a general initial 
condition, by linearity. However, such a superposition can only be composed for 
the non-centralized moments. At large times the first and second moment appear 
to increase at a constant rate, which suggest that there exist an effective velocity 
Ve and effective dispersion coefficient Df,: 

A + ^ (A + ^j'^ A + /_* 
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5.2. Two-dimensional. The 2D case is described by the equations (16). We may 
distinguish two different type of moments. The first type is that of the marginal 
moments, or moments with respect to the a;-coordinate ignoring the information on 
the ^-coordinates of the particles. The second type consist of conditional moments, 
either a moment with respect to x for a given value of y, or a moment with respect 
to the y and for a given value of x. Interestingly, the marginal moments for the 2D 
case are identical to the moments for the ID case. In the following subsections we 
discuss the conditional moments. 

5.2.1. The x-moments conditioned on y. This category of moments represent ex- 
pected values of x" for a population of particles with a specific y-coordinate. Note 
that these moments are a function of y. We define these moments as: 

-t-oo +00 

Mf{y)^ J x^Nf{x,y)dx, and Aff ) (y) = J x^Na{x,y)dx 
—00 —00 
For n = the expressions represent for each phase the particle distribution along 
the y-direction (the total mass of particles with the specified y-coordinate) . Note 
that the higher order moments {n > 0) are not yet divided by the zeroth moment, 
so higher order moments are not normalized. For the initial condition we consider 
a (unit) pulse with the phases in equilibrium. Then, the two-dimensional particle 
distribution becomes: 



where i is a or /. This expression is obtained by applying (20 1, replace hij by 



/i^* (see (15)) and take equal to 1. For the equilibrium initial condition the 
moments are: 

-t-00 



M^iy)^ / x^Nl\x,y)dx 



When we use (21) and change the order of integration, the zeroth, first and second 



moments become: 



(22) Mrfa)=/*/.:'(.,.)^d. 

(23) M'>'\y) ^ [ vrhrir,t),^^^dr 

Jo 2^/'kDt t 

(24) Mf)(y)=^(2Dir+t;V2)/if(r,t)^^dr. 

Figure [6] shows the distribution along the y-axis for the zeroth moments of free 
and adsorbed phase for several times. The initial condition here is the equilibrium 
situation. The zeroth moment is plotted horizontally. The figure shows that par- 
ticles of both phases gradually spread out in y-direction. Note in the figure at the 
right hand side that at t—^ the adsorbed phase is concentrated along y = and 
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Figure 6. Zeroth moments as given by (22 1 versus y for t=5,20,80, 
left free phase, right adsorbed phase; X — ^ — 0.05; Dt =0.03. 
Zeroth moment on the horizontal axis. In the figure at the right 
hand side the curve for t=5 has a pulse-shaped component 5{y) at 
y=0. 



the curve has a pulse-shaped component 5{y) for ?;=0. This pulse represents the 
particles that did not yet spend time in the free phase and are still in the initial 
position. At larger times most adsorbed particles do have spent time in the free 
phase and the spreading in y-direction becomes visible. 

Figure [7] shows the first moment for the free phase. The first moment, plotted 
horizontally, can be interpreted as the average distance travelled in the x-direction 
for the particles with a specific y-coordinate. As it appears, for ?/ the first mo- 
ments are smaller than for values of y greater or smaller than 0. It is an alternative 
illustration of the tailing effect. Particles spending less time in the free phase have 
less opportunity to displace in the x-direction and spread in the y-direction. Thus, 
a considerable fraction of these particles are found around the point of origin. Fur- 
ther away from the y-axis particles occur that have been able to disperse laterally. 
Therefore, they did spend some time in the free phase and, consequently, were also 
displaced further along the x-axis. Once more, we conclude that displacements 
in X- and y-direction are mutually dependent. This dependence in the spreading 
pattern is non-Gaussian. 

5.2.2. The y-moments, conditioned on x. The y-moments conditioned on x are 
defined as: 

oo oo 

M^p\x)^ j y^Nf{x,y)dy, and Mi^\x) ^ j y^N,{x,y) Ay, 

— oo —oo 

We analyze the y-moments by the Laplace Transform using the following initial 
conditions: 

Nf{x,y,t)\^^Q^TTf5{x)5{y), and Na{x,y,t)\^^Q ^ TTaS{x)6{y) 
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Figure 7. First moments for free phase versus y from (23 1 for 
t = 5, 10, 20. Numerical values: A = 0.05; n = 0.05; v = 0.3; 
Dl — 0.3; Dt = 0.03. First moment on the horizontal axis. 



After applying the Laplace transform to eq (161 and eliminating the transform of 
Na we obtain the following differential equation: 

+ - - '^^f + b.f5i.)5{y) = 

where Nf denotes the particle density in Laplace space, s is the Laplace parameter 
and b is an expression depending on s: 

s + X + fj. 

b = 

s + fi 



Na is related to Nf by (see eq. ([l6|)): (s + iJ.)Na - TTaS{x)S{y) = XNf. 

After taking the moments with respect to y we find for the zeroth moment: 



The solution for My' = My' (x) is 



^^^^ - " -'^W + W^(^) - 

A 

s + fi s + fl 

(0) _ iCFio), 



Miy» ^ -^Mf + ^5{x) 



jg(0) ^ hiTf / xv-\xyv^ +AbsDL 

where we use that the concentrations are zero at a; = ±oo. The inverse can be found 
with numerical Laplace inversion algorithms [5] , [32] , [I] , and is shown in Figure [s] 
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Figure 8. Zeroth moment for total concentration M^") (dashed), 
free concentration il/j"-* (dotted) and adsorbed Aia^^ concentration 
(drawn). Small numerical inaccuracies lead to the fact that the 
zeroth moment for the total concentration is below the curve for the 
free concentration. The peaks at a; = 0, representing the function 
S (x) are "approximated" by two exponential curves for reasons of 
illustration. 



For the second moment we obtain the equations: 

Dl—4 v^:-^ sbMf^ + 2DtM^°^ = 0, and sM^ ^ aM^^ - ^iMP 

OX'' ox ^ ^ ' 

oo 

where we use that j ^~g^ = Afj"'' , with the solution: 

— oo 

^(2), , , fl^DL'^'^Py 2DI ) 2Di 



MY'ix) = bTTf^^^ ^— -f^ x^{x+ , " , for X > 0, 

and Mi^\x) = -^Mf\x). 

s + 11 ' 

The ratio of the numerically [32] inverted m'^P and Mj.^-* is the normalized and 
centralized second moment (conditioned on x) and this ratio turns out to be x/Pe. 
Thus, the lateral spread depends on x, which (again) explains the tailing and the 
previously observed 'cigar shape'. 



6. Conclusions 

We discuss a stochastic particle approach, considering a particle that changes 
between a mobile and immobile state governed by a Markov chain, while its spatial 
displacement is governed by a random walk. Generally, particle models are assumed 
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to describe the advective-dispersive-kinetic transport process correctly, but in the 
hterature a rigorous proof is lacking. Our analysis shows that it is possible to derive 
the correct set of differential equations from a stochastic model for a single particle. 
To examine the non-Gaussian nature of the spreading process we analyze first the 
telegrapher's equation. This equation arises when only advection and kinetic sorp- 
tion is considered. We show that in such a system an apparent dispersion process 
occurs, generated only by the kinetic changing of the particle in states with different 
velocities (i.e. zero or v). This 'kinetics- induced' dispersion is non-Gaussian for 
short and intermediate times, while at large times the process develops as Gaussian 
dispersion. When hydrodynamic dispersion is included again, the spreading process 
becomes a combination of a Gaussian and non-Gaussian dispersion. We illustrate 
this for the 2D case. At short times the process is Gaussian, since hydrodynamic 
dispersion is the dominating process. At intermediate times the influence of kinetic 
induced dispersion increases and the spreading becomes non-Gaussian. Finally, at 
large times the dispersion becomes Gaussian again, but the (effective) longitudi- 
nal dispersion coefficient has an additional term due to the kinetics. Moveover, 
we find for the 2D case that the transverse spreading depends on the longitudinal 
coordinate, resulting in 'cigar-shaped' contours. The mechanism is best illustrated 
when longitudinal dispersion is assumed zero. Here the particles displace in the 
x-direction by advection and spread transversely by dispersion. Particles spending 
more time in the adsorbed phase are displaced less in x-direction, and also less 
spread out in y-direction. In a truly Gaussian distribution the transverse spread- 
ing is independent of the longitudinal coordinate. When longitudinal dispersion 
is included the same effect is observed, although for short times (compared to the 
kinetic exchange rate) the situation is now dominated by hydrodynamic dispersion. 
With respect of the validity of effective properties (velocity and dispersion) , we con- 
clude the following. The velocity and dispersion coefficients are represented by the 
rate of increase of the first and centralized second moment (times 1/2). For cases 
with low adsorption and desorption rates, the rates of increase for the moments re- 
main time dependent for a relatively long time. We conclude that constant effective 
properties can not be defined directly after the start of solute injection. For large 
times an asymptotic behavior is observed with a constant mean displacement and 
rate of dispersion, while the third moment vanishes and the kurtosis approaches a 
value of 3. This can be proved via our stochastic model by applying a sophisticated 
version of the Gentral Limit Theorem (Section 2.2). The critical time, required 
before dispersion coefficients becomes constant, is in the order of 3/(A -I- ^). The 
effective velocity is vfi/{X + similar to the case of linear equilibrium adsorption. 
However, at early times the free particles move with the original groundwater ve- 
locity and build up a lead with respect to the adsorbed phase. In the asymptotic 
stage, the adsorbed and free particles displace with the same average velocity, but 
the lead of the free plume is maintained. The effective longitudinal dispersion, can 
be much higher than in the case of linear equilibrium adsorption. The additional 
term, t;^A/i/(A + /i)'^ depends also on the groundwater velocity. It is remarkable 
that it depends on u^, while the micro-scale dispersion is linear in v. 
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